#import libraries
library(readr)
library(dplyr)
library(tidyr)
library(tidyverse)
library(ordinal)
library(effectsize)
library(emmeans)


#load data
data <- read_csv("Study 4b.csv", show_col_types = FALSE)


#recode important categorical vars with values
data$gender_coded <- recode_factor(data$gender, 
                                   '0' = 'male', 
                                   '1' = 'female',
                                   '2' = 'non-binary or other',
                                   .default = 'missing')
colnames(data)


#get the sample information
##condition
condition_count <- table(data$condition) #high_inequality_high_ses = 80,high_inequality_low_ses = 80, high_inequality_no_ses = 80, low_inequality_no_ses = 80
condition_count 

##gender
gender_count <- table(data$gender_coded) #male = 148, female = 169, non-binary or other = 3
gender_count 

##age
mean(data$age) #M = 42.98
sd(data$age) #SD = 11.19


#ANOVA
##anova
model <- aov(willingness ~ condition, data = data) #F(3, 316) = 18.11 (p < .001)
summary(model)

##get effect size
effectsize:: eta_squared(model, partial = TRUE) #eta2 = .15

##get descriptive
###get raw means
aggregate(willingness ~ condition, data = data, mean) #high_inequality_high_ses = 5.66, high_inequality_low_ses = 4.45, high_inequality_no_ses = 5.19, low_inequality_no_ses = 6.03
aggregate(willingness ~ condition, data = data, sd) #high_inequality_high_ses = 1.29, high_inequality_low_ses = 1.73, high_inequality_no_ses = 1.51, low_inequality_no_ses = 1.10

###get marginal means & contrasts
emm <- (emmeans(model, ~ condition))
pairs(emm, adjust = 'None')
eff_size(emm, sigma = sigma(model), edf = df.residual(model))